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Abstract 

Traditional epidemic detection algorithms make decisions using only local 
information. We propose a novel approach that explicitly models spatial 
information fusion from several metapopulations. Our method also takes 
into account cost-beneht considerations regarding the announcement of epi¬ 
demic. We utilize a compartmental stochastic model within a Bayesian 
detection framework which leads to a dynamic optimization problem. The 
resulting adaptive, non-parametric detection strategy optimally balances de¬ 
tection delay vis-a-vis probability of false alarms. Taking advantage of the 
underlying state-space structure, we represent the stopping rule in terms of 
a detection map which visualizes the relationship between the multivariate 
system state and policy making. It also allows us to obtain an efficient 
simulation-based solution algorithm that is based on the Sequential Regres¬ 
sion Monte Carlo (SRMC) approach of Gramacy and Ludkovski (SIFIN, 
2015). We illustrate our results on synthetic examples and also quantify 
the advantages of our adaptive detection relative to conventional threshold- 
based strategies. 

Keywords: Biosurveillance; quickest detection; regression Monte Carlo; 
stochastic compartmental models; 


1. Introduction 

Infectious disease epidemics intrinsically unfold across both space and 
time. As a result, biosurveillance algorithms need to integrate spatio-temporal 
data. This is especially so in the context of statistical inference, whereby 


‘Corresponding author 

Email addresses: shatskikh@pstat.ucsb.edu (Katherine Shatskikh), 
ludkovski@pstat.ucsb.edu (Michael Ludkovski) 


Preprint submitted to Elsevier 


September 15, 2015 



syndromic surveillance at neighboring locales carries additional information 
that can be fused for improved decision making in terms of initiating and or¬ 
ganizing epidemic counter-measures. A crucial first step for response strate¬ 
gies is to identify, or detect, in real-time the epidemic outset. In this article, 
we propose a methodology that allows for such optimal decision-making 
with spatial information fusion. Specifically, we investigate a model that 
combines quickest detection with a spatial metapopulation setup, integrat¬ 
ing information received from multiple geographic domains. To reflect the 
inherent uncertainty in epidemic evolution (which is amplified under partial 
information), we develop a stochastic compartmental (or state-space) epi¬ 
demic model, which allows us to generate adaptive, nonparametric detec¬ 
tion rules. Extant approaches largely propose heuristic detection strategies, 
concentrating primarily on the inferential aspect of the statistical model 
[HEIEIII!. For instance, a typical approach is to announce an epidemic as 
soon as the estimated number of infecteds in the local population is above a 
fixed I. In contrast, we dynamically optimize the detection strategy, to come 
up with a “best” detection rule within our mechanistic outbreak model. 

Traditional compartmental epidemic models deal with a single popu¬ 
lation; the spatial aspect is treated by building a series of such single¬ 
population models that are estimated/forecasted independently. This is 
also a common surveillance approach, especially for recurring infectious epi¬ 
demics, such as influenza-like illness (ILI), dengue fever, or measles. For 
example, in the US the existing biosurveillance systems for flu operate pri¬ 
marily at the state level and are siloed across states. This limitation of 
existing practice was brought into sharp relief during the 2014 Ebola out¬ 
break in West Africa. The epidemic has been accompanied by a dearth of 
reliable information, leading to extreme spread in forecasts regarding the 
future course of the outbreak. In addition, numerous statistical methods 
laEiiT! were put forth attempting to infer in “real-time” the actual size and 
parameters of the outbreak in different locales. However, nearly all these 
methods were single-population, so that when trying for example to infer the 
number of Ebola infecteds in Liberia, only Liberian data was utilized, com¬ 
pletely ignoring similar and highly relevant data from neighboring Guinea 
and Sierra Leone. Similarly, at the more granular provincial level, data from 
neighboring provinces was generally not used during estimation procedures. 

For a less dramatic and perhaps more statistically convenient example, 
we discuss the yearly influenza outbreaks in United States. Figure illus¬ 
trates the spatial dynamics of ILI during the 2012-13 flu season. As can be 
observed, the peak of the outbreak varied significantly (up to 6-8 weeks dif¬ 
ference) across different parts of the country. Nevertheless, there is a clear 
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propagation, making spatial information fusion desirable. Figure [^indicates 
that the current, single-population based detection protocols are not suffi¬ 
cient; for instance the fact that there are increased ILI levels in Arizona is 
ought to be taken into account when trying to detect or forecast the epi¬ 
demic start in California. A further important remark is that the illustrated 
spatial spread is year-specific, and in other years rather different patterns 
may be observed. 




Week 49, 2012 


Week 1, 2013 




Week 4, 2013 


Figure 1: Spread of Influenza during the 2012-13 Flu season according to 
FluView CDC data. The colors represent weekly ILI activity levels in terms 
of percentage of doctor visits attributed to ILI relative to low-season base¬ 
line. Green indicates at/below mean, while shades of red indicate outbreak 
activity (with darkest color corresponding to eight or more standard de¬ 
viations above the mean). Weeks are numbered from January 1, and are 
12/3-9/2012 (Week 49), 12/31/2012-1/6/2013 (Week 1) and 1/21-27/2013 
(Week 4), respectively. Data source: http://www.cdc.gov/flu/weekly/ 
pastreports.htm. 


1.1. Contributions 

In this paper we formulate and analyze an epidemic detection problem 
within a multi-population paradigm. To do so, we develop a reduced com- 
partmental model that extends the classical Susceptible-Infected-Recovered 
(SIR) setup to two population pools. Pools are interpreted as distinct geo¬ 
graphic regions, e.g. states or counties. To fix ideas, we consider the situation 
where the epidemic begins in Pool 1 and subsequently may be transmitted 
to Pool 2 via infecteds that travel between the two pools. The aim of the 
policy-maker is to detect, as soon as possible and in online fashion, the onset 
of epidemic in Pool 2. 

To capture the inferential aspect, we assume that full information is 
available about the outbreak in Pool 1, but only partial information about 
Pool 2. As a result, one has to make imperfect decisions and in particular 
address the canonical trade-off between making announcements too early 
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(so called “false alarms”) and making decisions too late (“detection delay”). 
Indeed, if the detection is too late, then a certain number of infections would 
be missed and it would be harder to stop the epidemic from spreading. If the 
detection is premature, human, financial and reputational resources would 
be wasted. Therefore, a careful trade-off between those costs should be done 
to balance costs due to epidemic morbidity and costs arising from policy 
actions. We then use the above cost analysis to quantify decision-making 
quality and to dehne optimality of detection strategies. 

Mathematically, we cast the online detection problem as a dynamic opti¬ 
mization problem, connecting to the classical dynamic programming formu¬ 
lation [S] in control theory. A major challenge with dynamic programming 
(which is perhaps the prime reason for the lack in its uptake in the bio¬ 
surveillance community) is computational bottlenecks due to the curse of 
dimensionality. Indeed, the above optimization problem is nontrivial from 
several directions. First, because the underlying system is stochastic, the 
optimal solution is adaptive, i.e. a function of the current system state. 
Consequently, there is no simple description to the resulting detection strat¬ 
egy which is instead summarized through a detection map that translates 
system states into optimal detection decisions. Second, the nonlinear dy¬ 
namics of the SIR model preclude analytic solutions. Crucially, there are no 
analytic expressions for the future distribution of the system state, which 
necessitates the use of numerical approximations to solve the optimization 
problem. Third, because the system state is multivariate and too large to 
enumerate, the corresponding integrals are computationally demanding. 

However, taking advantage of the detection strategy structure, which 
requires simply announcing at each stage whether the epidemic has reached 
Pool 2 or not, we implement an efficient numerical algorithm. Specihcally, 
we rely on the recent Sequential Regression Monte Carlo (SRMC) method 
of [9j, which blends modern statistical tools, including nonparametric re¬ 
gression and sequential design, with approximate dynamic programming, to 
drastically mitigate issues of computational efficiency. 

The main contributions of this work are then threefold. First, we propose 
and analyze a multi-population extension of the classical SIR model, as 
well as a reduced version suitable for the Bayesian detection framework. 
Second, we develop and adapt an extension of the sequential regression 
Monte Carlo (SRMC) approach to efficiently solve the dynamic optimization 
problem. Third, we present a detailed investigation into the performance of 
the designed strategy, in particular in comparison to conventional threshold- 
based strategies. 

The organization of the paper is as follows. Section formalizes the 
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mathematical aspects of our model, including the detection setup. The 
stochastic dynamics of the outbreak are rigorized in Section Section 
presents numerical illustrations of our method as well as comparison with 
other methods. Section then describes the Sequential Regression Monte 
Carlo algorithm that we developed for our setup. Section provides the 
conclusion and the discussion on future extensions of our framework. 

1.2. Spatial Stochastic Epidemic Models 

Mathematical models of infectious disease epidemics have become an im¬ 
portant tool in the arsenal of public health policy. In an idealized world, 
detection reduces to the mathematical problem of clustering, tracking the 
health status of the surveyed individuals and identifying unusual aberra¬ 
tions in either the temporal or spatial dimensions. In reality, there is the 
additional aspect of missing information which necessitates the application 
of statistical inference algorithms, as well as a mathematical model for the 
epidemic. In the context of online inference, a simple mechanistic approach 
that allows for maximum tractability continues to be the most popular, and 
is also adopted here. Specihcally, we rely on the formalism of an SIR model 
|10j that implies proportional homogenous mixing between infecteds and 
susceptibles within a population pool. Spatial heterogeneity is captured by 
incorporating meta-populations, also known as patch models [HI [m dg. 
The multi-patch approach partitions the global population into distinct dis¬ 
crete regions or pools, allowing for local spread of the epidemic within each 
pool, as well as global transmission that is specified via a mobility matrix. 
As in mm we assume that susceptibles are stationary, while infecteds can 
move or travel between the pools, creating cross-infections. 

Alternative frameworks for epidemic spread include point process models 
El, and network models m that provide more nuanced interaction between 
individuals to mimic existing social structures, such as households, schools, 
and workplaces. At even more detail, agent-based models m generate 
micro-simulations that provide a detailed synthetic view for each individual 
and their social interactions. Such models can also incorporate precise travel 
patterns HH However, the latter paradigms are geared towards realistic 
forecasting of epidemic progress and are less suited for online detection due 
to intractable inference in terms of observed data and the computational 
expenses in generating micro-scenarios. 

A variety of approaches exist for constructing outbreak detection rules, 
see for example the recent survey by Shmueli and Burkom |T8j, and the 
monograph by Lawson m- Quality control methods m introduced in the 
1950s form the simplest class of rules and continue to be common. Other 
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heuristics include moving-average tools [2T], various scan statistics Eam, 
and branching-process approximations [53] . More explicit cost-benefit anal¬ 
ysis for the trade-off between false alarms and detection delay can be applied 
using the Cumulative Sum (CUSUM) framework |24] . CUSUM also under¬ 
lies the early aberration response system (EARS) employed by the Centers 
for Disease Control [25]. Alternatively, Bayesian methods allow to further 
assess the uncertainty involved in decision-making based on partial infor¬ 
mation. Two main types are hidden Markov models |M1I27| and Bayesian 
hierarchical models PEHI. The Bayesian paradigm translates epidemic data 
into the posterior probability of an outbreak. To convert the latter into a 
detection rule, one typically employs a simple threshold strategy. For ex¬ 
ample, in [T], the authors recommend “an alert for action if the posterior 
probability is larger than 70%”. We further refine this approach by de¬ 
riving optimal, non-parametric detection strategies based on the inputted 
cost-benefit parameters. 

Detection can be seen as a basic form of epidemic response, and indeed 
our computational methodology can be extended to this more general prob¬ 
lem. In that sense, this paper extends the first author’s previous work on 
stochastic control methods for controlling epidemics |29l ISO] . Similar to |29| , 
we design a Bayesian dynamic optimization algorithm for biosurveillance de¬ 
cision policy. Other mathematically oriented studies that consider optimal 
control of epidemics include |3I1E2|. 

In the context of detection with limited information, a spatial epidemic 
model requires information fusion. Fusion of information channels for the 
purpose of biosurveillance has been an area of intense research in the past 
decade. On the one hand, novel information sources, such as social me¬ 
dia |33| or internet data |3l| have created new opportunities for syndromic 
surveillance. On the other hand, developments in statistical fusion tech¬ 
niques [la ESI ESI have led to new ways of integrating multivariate infor¬ 
mation streams. In particular, there has been a lot of interest in online 
Bayesian approaches plElCSlEilEll that allow for predictive modeling and 
forecasting of epidemics. The above models all focus on a single homoge¬ 
nous population with the different surveillance channels complementing each 
other. In contrast, we consider multiple underlying population pools each 
with a distinct, but co-dependent information channel. In terms of ex¬ 
plicitly accounting for spatial propagation, our work is closest to [38] who 
considered a spatial “wave” model for an epidemic. In the present article, 
we connect this framework to the SIR context, modeling epidemic spread 
across geographically-based population pools. The resulting decision strat¬ 
egy provides insights into integrating data from multiple spatial locales for 
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the purposes of detection, cf. Section below. 

2. Quickest Detection 
2.1. Mathematical Model 

We work with a state-space model, denoting by Xt the epidemic state at 
times t = 0,l,2,.... A typical length of one time period in biosurveillance 
is a week. The precise components of X will be specihed later; abstractly 
X is taken to be a stochastic Markov process taking values in a state space 
X C M'^, and summarizes information about both Pool 1 and Pool 2. In 
particular, X contains information about the number of infecteds in 
Pool A: = 1,2 at time t. The transition kernel of X is assumed to be time¬ 
stationary and is denoted by ps(x|y) = P{Xt+s = x|Xi = y), x, y G X. 

The aim of the policy maker is to detect the onset of epidemic in Pool 2. 
A detection strategy is probabilistically represented as a dynamic “alarm” 
which announces an outbreak in Pool 2, based on information gathered so 
far. Only a single announcement is allowed; once announced, the detection 
problem is assumed to be over. The set of such detection strategies is ex¬ 
pressed through the set S of X-stopping times, where Xt = cr(Xo:t) is the 
information filtration generated by X by time t. A strategy r G 5 is a ran¬ 
dom variable taking values in r G {0,1,2,...}, such that {r = t} G Xt (this 
requirement captures the fact that r must be “online” in terms of the infor¬ 
mation available so far). Thanks to the Markov property of X, the structure 
of r can be summarized via a detection map. Indeed, at each time-step 
there is the binary decision to either “announce” an outbreak (subset 6), 
or wait for another period (subset Gl). Since the evolution of X is stationary 
in time, the corresponding partition of the state space is also independent 
of t. Dynamically, this implies that r announces the epidemic the first time 
that the state X enters the region © C X, 

T = inf{t : Xt G S}. (1) 

Equation Q gives a one-to-one correspondence between detection strategies 
r and detection maps &. In other words, the detection strategies we consider 
are of online feedback type, based on the trajectory of X. 

As mentioned, the dynamic optimization objective consists in optimally 
trading off the concern of premature announcements against any potential 
delays. These conflicting costs are measured through the immediate stop¬ 
ping cost d(x) and the cost of waiting. The immediate costs are linked to 
the penalty for false alarms, specified by a given constant Cfa- We assume 
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that CpA is paid if and only if the epidemic has not yet reached Pool 2, so 
that 


(i(xo) CpA • (2) 

Waiting costs are assumed to be proportional to detection delay, i.e. the 
time between the outbreak reaching Pool 2 and outbreak announcement. 
Define 6 to be the time when the second population gets infected from the 
first population, i.e. 

9 := inf{t : = 0 and > 0}. 

Then the detection delay is max(r—0,0) and carries cost Cpeiay max(r—0,0). 
This is equivalent to charging waiting costs of CDeiayl|j{ 2 )^Q| at each step 
until surveillance is terminated at the random instant r, so that total waiting 
costs on [0, r] are 


r—1 

c(Xo:r) := 'y ^ C'DeIayl|r(2)^p,| “I" ^Fa1|j{2)_q|. (3) 

s=0 

We will refer to the costs d{-) and c(-) as the immediate cost and the future 
cost, respectively. 

Remark 1. Note that detection costs are intrinsically defined in terms of the 
count of infecteds in Pool 2, which is assumed to be unavailable to the 
policy-maker. Below we will operationalize Q and Q by taking conditional 
expectation with respect to information that is available, see 

The aim of outbreak detection is to pinpoint 0, i.e. ideally one takes 
T = 6. However, this is not possible if only partial information is available 
about X, specifically about B . When r and 0 are different, Coeiay penalizes 
the event {r > 0}, and CpA penalizes {r < 0}. The cost structure in ([^ 
is then a dynamic counterpart of the usual Type-I and Type-II errors in 
hypothesis testing. 

2.2. Detection Problem 

Our detection problem is formalized as minimizing the expected future 
cost over all possible stopping times r [39], i.e. an optimal stopping problem. 
Namely, we define the value function V as 

V (xo) := inf E [c(Xo:r)|Xo = xq] , (4) 
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where xq is the initial state. Assuming the infimum in Q is achieved, the 
dynamic programming principle [SU] implies 

l/(xo) = min(d(xo),.E [l/(Xi)|Xo = xq]) , (5) 

where the conditional expectation operator is 

E[V{Xi)\Xo = xo] = y y(x)pi(x|xo)dx. 

The minimum operator in corresponds to the idea that it is optimal to 
declare an outbreak if the immediate cost is smaller than the future cost, 
i.e. the likelihood of false alarms is dominated by the cost of waiting. The 
former case is equivalent to the expectation of the value function at time 1 
being greater than the immediate cost, and therefore we may classify the 
stopping region via 

6:={x:£;[l/(Xi)|To = x]-d(x) >0}. (6) 

Hence, in terms of the above detection map, our goal is to optimally partition 
A = ©UGl into two regions, such that © consists of all initial states xq where 
it is optimal to declare the epidemic, and C is its complement, where it is 
optimal to wait. 


2.3. Reduction to a Model Predictive Control Problem 

The characterization in ^ is implicit, since it features H(-) on both sides 
of the expression. Specifically, the value function V corresponds to a fixed 
point [8] of the functional operator C, defined by (£u)(x) := min (d(x), £' [u(Xi)| To = x]). 
To solve for V (x), a basic strategy is then to apply Picard-type fixed-point 
iterations. In other words, given some initial guess p(^^(x), we build a 
sequence of approximations via := or explicitly. 


fW(xo) 


mm I 


(d(xo),P;[H('=-^)(Ti)|To 



( 7 ) 


However, to guarantee the convergence of does not appear tractable, 
and the practical performance of Q is very sensitive to the initial guess 
To circumvent this challenge, we rely on the concept of model predic¬ 
tive control. To wit, we introduce an auxiliary parameter t which can be 
intuitively thought of as forward time. The value functions V{t,-) and de¬ 
tection maps &t are now also indexed by t. We start with the trivial initial 
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condition 1/(0,x) := d{x), which corresponds to ©o = Next, mimicking 


the classical dynamic programming on finite horizon, we define 

l/(f,xo) := min(d(xo),-E [y(f - l,Xi)|Xo = xo]), t = l,2,.... (8) 

Define the Q-value, also known as costs-to-go by 

q{t,x):=E[V{t-l,Xi)\Xo=x]. (9) 

Then the stopping set at iteration t is 

©t := {xo G Af : g(f,xo) - d(xo) > 0} , t = l,2,.... (10) 

We may “unroll” the expectation encoded mV{t — l,Xi) to write 

q{t,xo) = E[V{t- l,Xi)|Xo = xo] = [c(Xo,^(t))|To = xq] , (11) 


where = min{s > 1 : G ©^-s}. This justifies the interpretation of 

q{t,-) as costs-to-go, since c(Xq.^{o) are indeed the future costs associated 
with not stopping immediately. 



Figure 2: Detection strategy at iteration t = 1. The example is based on 
the model of Section]^ with parameters in TableIn the plot, the state of 
Pool 1 is held fixed at 5*0^^ = 1990, = 10. 


Figure illustrates the first step of the recursion Q at f = 1. In the 
plot we compare 


^[F(0,Xi)|Xo = x] =^[d(Xi)|Xo = x] 
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against d{x). As discussed, the epidemic is announced when E[V{t — 1, j£i)|Xo = x] > 
d{x) (the right side of the plot). In the opposite case (the left side of the 
plot), the optimal decision is to wait. As shown by the Figure, the structure 
of the decision map is driven by the regions where these two quantities are 
equal to each other, which corresponds to the detection boundary, 


d&t := {x : S [V{t - 1, Xi)|Xo = x] = d(x)} 


( 12 ) 


The stopping region &t is our detection rule at iteration step t. It can 
be characterized as the optimal detection rule among all strategies in 5^*^ = 
{r G T" : r < t} that are upper-bounded by t (By construction, < t). 
As t —)> oo, we have that the set of admissible rules expands 5^*^ S, and 
hence we expect that ©t —)■ S and V{t,x) —)• ^(x) . Intuitively, for large t, 
the recursively defined Q converges to a stationary case that ought to be 
the fixed point defining V (x) in Q . The above convergence can be improved 
via model predictive control (also known as receding horizon control) [40j 
which applies the fixed dete ctio n map rather than the time-dependent 
6t at each step, cf. Section 5.1 


3. Epidemic Model 

3.1. Multiple Population SIR model 

A susceptible-infected-recovered (SIR) model provides an aggregate “grav¬ 
ity” view of the epidemic by focusing on three basic types of individuals in 
the population: susceptible, infected and recovered. Susceptible individuals 
are the ones who haven’t experienced the disease yet. Interaction between 
an infected and susceptible individuals can lead to an infection. Thus, con¬ 
tacts stochastically generate new infecteds who in turn can further infect 
other susceptible individuals. After some time an infected individual recov¬ 
ers and becomes immune (i.e. becomes a Recovered): he/she can no longer 
infect others or get infected. 

While the detection problem is specified at the discrete instances t = 
1 , 2 ,..., for describing outbreak dynamics it is more convenient to work with 
continuous-time dynamical systems. As in m Ch. 6], we thus first recall the 
multi-type stochastic SIR model in continuous time. The overall epidemic 
state at epoch t G ]R_|_ is denoted by {S*, R, R* }, where St = ..., 

It = ..., and Rt = {r[^\ ..., are vectors denoting the 

count of susceptible, infected and recovered individuals in each of 1 < A: < A' 
meta-populations. We assume that the pool size of each meta-population is 
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fixed at As a result, we omit further mention of 

since it can be found from — l[^\ 

The continuous evolution of the state process {St, It} G {(s, i) : + ifc < 

Vfc < K} is described through Markov chain or stochastic kinetic 
system language. Namely, the epidemic state is piecewise constant in time. 
Next, there are 2ii'+A'(iir —1) possible transitions, described by the reaction 
channels |41j : 

Infection 5(*^) +/(*=) w/rate 

< Transmission /(^l _)_/(^') w/rate * 

Recovery —)• 0 w/rate 

(13) 

where each reaction is further indexed by 1 < A: < AT and k' < K, k' ^ k. 

The first transition represents an infection of a susceptible individual by 
an infected individual from the same pool k. This transition happens at 

(k) 

rate where 1 < fc < AT and I3k is a contact rate of infected and 

susceptible individuals within the k-th meta-population. 

The second transition is a transmission: an infection of a susceptible 
individual from pool k by an infected individual from a different pool k'. 

(k') / 

The frequency of such infections is I3k,k'lt where 1 < k,k < K and 

(3k^k' is a contact rate of infected and susceptible individuals from different 
populations. Since contacts between individuals from different populations 
are less frequent, /3fc,fc' (3k'- To reduce the number of parameters, we thus 
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assume that cross-population interactions occur at rate I3k^k' = oth'-, where 
a is the proportion of “travelers” in each pool. Thus, cross-contacts happen 
at the fraction a of a contact rate within one population; a typical range is 
a G [0.01,0.2], 

The last transition in (13) is a recovery and subsequent immunity of an 

* - ' (fc) 

infected individual in a population k. The rate of transition is 7 /^ , where 7 

is a recovery rate, independent of the pool index k. This can be interpreted 

as individuals staying infected for an Exponentially distributed time with 

mean I/ 7 . 

In this paper we focus on two-population models, positing that the out¬ 
break begins in Pool 1 and may subsequently spread to Pool 2, where it 
is to be detected. Accordingly, we will be fusing information from Pool 1 
and Pool 2 to identify the onset of epidemic in Pool 2. We assume that the 
two pools have similar characteristics, so that all parameters are homoge¬ 
nous in A: = 1,2. Thus, the two-population SIR model, shown in Figure 
has 3 parameters a - the mixing parameter between two populations, /3 - 
the within-pool contact rate of infected and susceptible individuals, and 7 
- the recovery rate. These parameters are assumed to be known and are a 
function of the modeled disease family (e.g. influenza or dengue fever), the 
demographics and public health characteristics of the populations, and the 
travel patterns across pools. 


3.2. Partial Observation 

Under full observations the detection problem Q would be trivial, since 
one can directly track ' and declare an outbreak as soon as there any 
infecteds in the second pool. However, realistically is not observed. 
Some of the reasons include mis-diagnoses among infecteds, patients not 
seeking care, false positives, mis-reporting or lack of reporting of epidemio¬ 
logical data, etc. Consequently, we assume that the true size of the S/I/R 
compartments in Pool 2 is not known. To simplify the presentation, we as¬ 
sume that is observed in Pool 1, perhaps due to better epidemiological 
surveillance in that pool. 

In our detection problem, the main event of interest is the presence of 
any infecteds in Pool 2, {I^ > 0}. Accordingly, we consider Pt = P{It > 

Oj^t), the posterior probability that the epidemic started in the second pop¬ 
ulation given the limited knowledge about it available by time t, here sum¬ 
marized by some information set Qt. Depending on assumptions about the 
observations structure, Pt may be available in closed form (e.g. through 
Bayesian conjugate updating [30l[2]) or may have to be only approximately 
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computed through e.g. particle filtering methods The latter method, 

/o'! 

which computes the whole posterior distribution vr^ ~ \Gt, is computa¬ 
tionally expensive, while conjugate updating requires carrying several suf- 
ficient statistics about the posterior of . In either case, Pt on its own 
is not Markovian, and hence does not possess simple dynamics. Therefore 
we propose a model that works with a simplified, Markovian version of Pt, 
which we denote as Pt- 



t, Time 


Figure 4: Posterior probability that the epidemic started in the second pop¬ 
ulation Pt vs. time t, where the dashed line represents the actual start time 9 
of outbreak in Pool 2. The plot was constructed using particle filtering using 
the following model parameter values: /3 = 0.75, a = 0.01, 7 = 0.5, M^-) = 
M(2) = 2000. 


Figure shows a sample scenario of the evolution of Pt in a partially 
observed framework. The plot was generated using particle filtering and 
used the two-pool model (13) with noisy Poisson-type observations in each 
pool 


We observe that Pt tends to drift up (i.e. posterior probability of 
outbreak increases over time) and eventually hits 1 . 


3.3. Reduced Model 

Our reduced model consists of the state of epidemic in the first popula¬ 
tion and a process Pt that is interpreted as the probability that 

the epidemic reached Pool 2 conditional on the information Gt = ^ ^o^t) 
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from Pool 1. The first two components and come from a one- 
population SIR model (see the definition of SIR model for R' = 1 in Sec¬ 
tion 3.1). To prescribe the dynamics of the pseudo-posterior we decom- 
pose the event {I^ > 0} = {0 < t} into two cases: the event that the 
epidemic already started at time t — 1 (i.e. 0 < t — 1), and the event that 
it starts at t = 9. We also add some stochastic noise to denote exogenous 
fluctuations in our posterior estimates regarding the second pool. In total, 
we thus assume that 


Pt = Pt-i + P(/S = 0 and 4"^ > 0\gt-i) + 6t, (14) 


where 5t are i.i.d. noise terms. Intuitively, the probability of outbreak has a 
positive drift over time, and the drift is precisely the posterior probability 
of the outbreak beginning during the current period, {6 £ [t — l,t]}. 

From the SIR dynamics, the probability that {0 G [t — l,t]} conditional 
on Pool-1 observations up to previous stage t — 1, is equal to the probability 
that an infected from Pool 1 interacts with a susceptible from Pool 2, times 
the conditional probability that {9 > t — 1}. The former happens with 
fl) 

rate a/3Is ■s £ — Ip], while the latter event is the complement of 

{9 < t—1} and hence has probability 1—Pt-i. Using the fact that conditional 
on {9 > t}, and making the transition rate constant on [t — 1, t] 

we obtain 

P{9 £[t- l,t]\gt-i) ~ a/3/S(l - Pt-i)- (15) 

To guarantee Pt £ [0,1] is interpretable as probability we confine it to 
[0,1], yielding 


0 V {Pt-i + - Pt-i) + 6t) A 1, 

1 , 


if Pt-i / 1 

if Pt_i = 1 


(16) 


In our simulations we use centered Gaussian noise 6t AA(0, cr|) with 
variance however it can take any distribution. Note that Pt = 1 is an 
absorbing state, representing certainty that the outbreak reached Pool 2, 
while Pt = 0 is a boundary case, since even if it is certain that the outbreak 
is currently not in Pool 2, it can still get cross-infected in the future. Similar 
features hold for the true posterior probability Pt, cf. FigureAlternative 
models for probability of outbreak Pt , are discussed in Section 


Remark 2. Note that (16) is in discrete-time; to connect to the continuous¬ 
time dynamics of SIR one could take the limit as the time increment goes to 
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zero, obtaining a diffusive model dPt = a/3ll:^\l — Pt) dt + ddWt where (Wt) 
is a Brownian motion. However, since detection is assumed to take place 
only at instances t = 1,2,..., we prefer to work with (16) as is. 


3 . 4 . Detection within the Reduced Model 

To sum up, the developed reduced 2-pool model has a 3-dimensional 
state {X}t = (sj:^\ ll:^\ Pt) with state space 


X := {{s,i,p) : s,i £ N, s + i < M^^\p £ [0,1]}. 

Figure shows a few sample trajectories of X to illustrate the resulting 
dynamics. 




t t 

Figure 5: Three sample trajectories of X with the initial condition = 
1995, = 5, Pq = 0 and outbreak parameters from Table Left panel is 

the plot of the number of infecteds in the first population, and right 

panel is the plot of {Pt}, the posterior probability that the epidemic started 
in the second population. The vertical dotted lines represent times when Pt 
hits 1 and outbreak becomes certain. 

Our detection problem ( |10| ) relies on the computation of the immediate 
and future expected costs P [c(jCo:r)|To] and d{Xo). Re-writing the def¬ 
initions of immediate and future costs ([^ and ([^ in terms of the event 
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( 17 ) 


( 2 ^\ 

{ir ^ O}? cLiid tciking conditionSil ©xpcctcition wg obtSiim 
d{Xo) :=Cfa(1-Po), 

r—1 

c(Xo;r) := ^ C'oelay-F’s + ^faCI “ Pt), 


(18) 


s=0 


where t € S. Rather than in terms of the unobserved the above 
expressions are now given in terms of the component Pt, allowing to measure 
detection costs within the X-model. Notice that d{Xo) is a function of Pq 
and c(Xo;r) is a function of the future trajectory {P^, s = 0 ,..., r}. 

Our goal is to find the detection maps ©t for t = 1, 2,..., defined recur¬ 


sively in (10). To do so, at each step we need to evaluate P [R(t — 1, Xi)|Xo = x] 
and fi(x). The immediate cost d{x.) can be computed exactly via How¬ 
ever, the expectation E[V{t — l,Xi)|Xo = x] can not be computed analyti¬ 
cally since there are no closed-form expressions for the distribution of Xo:r- 


In Section 5.3 we present the sequential Regression Monte Carlo approach 
which offers an efficient way to empirically estimate ©t based on syntheti¬ 
cally generated epidemic scenarios. We then use Model Predictive Control 
to estimate the stationary detection map ©. 


4. Case Study 

To illustrate the dynamic detection strategy within our 2-pool model, 
in this section we present a detailed case study. Table summarizes the 
parameters used. Epidemic parameters are taken to be /3 = 0.75 and 7 = 0.5. 
Thus, the initial reproduction ratio is TZo = /d/'j = 1.5, which is a moderately 
infectious epidemic. We assume that the pool mixing parameter is a = 0.01, 
which is reasonable for pools representing well-separated cities or counties. 


The inference noise in (16) is taken to be Gaussian with variance 6 t ~ 
= 0.01^). For the detection costs in ([I7l)-([l8l), we take without loss 
of generality Coeiay = 1 and fix Cfa = 20 . As we will see, this corresponds 
to a moderate penalty for false alarms. 


Epidemic: 

= 2000 

= M( 1 ) - 4^^ 

0-5 = 1/100 


/3 = 0.75 

a = 0.01 

7 = 0.5 

Costs/Penalties: 

> 

II 

to 

0 

C*Delay ~ 1 



Table 1: Outbreak and costs parameters for the case study of Section B crs 
refers to the noise in P, cf. (16). 
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Figure 6: Left panel: detection rule &20 ™ terms of and P. The 
detection boundary d& 2 Q is shown with the solid curve. We also show the 
experimental design Z that was used, illustrated with the scatterplot. Size of 
pixels corresponds to the number of times that neighborhood was sampled. 
Right panel: standard errors 'O(x) from (27). Observe lower standard errors 
in regions where the design Z is more dense. 


So far the case study features a three-dimensional state 
so that the resulting detection maps are in 3-D. To aid visualization, we 
consider a variant with a reduced dimension. Namely, we drop the com¬ 
ponent measuring the number of infecteds in Pool 1. Indeed, at the 
early stages of the outbreak the ratio is approximately one. As a 

result, one may assume that the rate of infections in Pool 1 is simply 
which corresponds to the classical branching process epidemic model [inj. It 
is known [l3] that this approximation remains valid up to t = 0(log(M(^^) 
by which time, = 0(VMW).; therefore it works especially well in large 
populations, and hence is termed a large-population (LP) approximation. 
The LP model only has two dimensions, X' := P} allowing to plot the 
corresponding 2-D stopping set 6^^. 

Figure shows generated under the conditions of Table and the 
above large population assumption. As expected, epidemic detection is trig- 
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gered once the posterior probability Pt of {/^ > 0}, is high enough. How¬ 

ever, we observe that detection is also highly sensitive to values of 
for instance detection is progressively delayed as gets bigger. This de¬ 
pendence between the two pools in terms of decision making illustrates the 
underlying cross-pool information fusion. Intuitively, detection should take 
place once Pt is high enough. However, conditional on a fixed Pt, larger 
number of Pool 1 infecteds makes an impending outbreak in Pool 2 more 
likely, lowering waiting costs. Hence, the detection boundary curves in 
Mathematically, recall that in (16), the growth rate of P increases in 
As a result, for large values of /AT one may expect that the next-stage Pt+i 
will also be large, i.e. move into the “Announce” region quicker. This again 
lowers the waiting costs and therefore delays announcement. 


4 . 1 . Evaluating Detection Rules 

Figure shows dynamic decision-making in the LP model through a 
collection of generated trajectories of X' = {lj:^\Pt} and their corresponding 
detection times the first time the state process X' enters the stopping 
set 6^^. We observe that the trajectories generally move north-east, as 
both P and tend to increase. However, the rate at which they grow and 
the precise direction are uncertain and vary across scenarios. Consequently, 
at detection, both P^lp and have a nontrivial distribution. 

To better understand the detection map , we analyze the resulting 
detection strategy given by and compare it to alternatives. Two classes 
of simpler detection rules are Threshold-P and Threshold-t. The Threshold- 
P strategy announces an outbreak as soon as Pt > P for a given threshold 
P. Hence, it acts solely based on local (posterior) information about Pool 
2. This mimics the CDC policy [53] of announcing an epidemic when the 
number of infecteds in Population 2 crosses some pre-specified level. In 
contrast to the fused detection strategy with a curved detection boundary 
which jointly takes into account both Pt and Threshold-P rule only uses 
Pt for detection decisions, yielding a flat, horizontal detection boundary in 
Figure The threshold-t strategy is a simple non-adaptive strategy that 
announces at the fixed stage t. It is illustrated in Figure where we record 
the joint distribution of , Pf at t = 8. 

Returning to the full 3-D model with state X we evaluate the result¬ 
ing optimal detection strategy t* and proceed to compare its performance 
against the other potential detection rules discussed above. Specifically, the 
first two alternatives are a Threshold-P rule with P = 0.8 (declare an epi¬ 
demic if its probability is above 80%) and a Threshold-t strategy with t = 8. 
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Number of Infecteds in Pop. 1 


Figure 7: Fifty sampled epidemic trajectories ema¬ 

nating from the initial state = 10 and Pq = 0.1. We show the LP 
detection boundary (namely dG^^), as well as a threshold strategy that 
announces epidemic as soon as P* > .P = 0.8. Lastly, the red crosses denote 
the locations of the trajectories at t = 8, which is the basis of the alternate 
Threshold-t strategy. 


The latter was found to be the best strategy among those that declare out¬ 
break at a fixed stage. The last alternative is the LP strategy from last 
section. Recall that makes decisions while ignoring In that sense, 
when applied to the full 3-D model, it gives a simplified, but still adaptive, 
detection rule. To recap, Threshold-t strategy is completely non-adaptive; 
Threshold-P only relies on P*; LP relies on Pt}, and Optimal strategy 
uses all of P^}. 

To compare the performance of the above competing strategies, we fixed 



Detection time r 
Mean StDev. 

Realized Cost Q 
Mean StDev. 

PFA P[1 - P^] 

Optimal 

8.86 

2.59 

6.53 

1.70 

8.2% 

LP 

9.32 

2.95 

6.57 

1.81 

6.4% 

Threshold-P 

7.88 

2.85 

7.03 

1.58 

15.3% 

Threshold-t 

8.00 

N/A 

7.18 

2.21 

14.4% 


Table 2: Comparison of Optimal, Large Population(LP), Threshold-P with 
P = 0.8 and Threshold-t with t = 8 strategies. Statistics are based on 1000 
synthetic trajectories of P}, where Q = c{XQ.^(t)). 
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Figure 8; Summary statistics of different detection strategies constructed 
from 1000 sample epidemic trajectories. The LP detection strategy is from 
Figure Right: Distribution of detection times r; Left: Distribution of 
posterior probability of outbreak in Pool 2 at detection time, Pr- 

the initial condition at S'q^^ = 1990, = 10 and Pq = O-lj so that there are 

10 infecteds in Pool 1 and 10% prior probability of epidemic already in Pool 
2. Then we simulated 1000 epidemic trajectories {xg.^}, n = 1,..., 1000, 
emanating from this fixed initial condition up to the detection time r (which 
depends in turn on the strategy used). Table then presents the resulting 
summary statistics based on these frozen 1000 trajectories (note that there 
are no analytic formulas to obtain these metrics, so we have to resort to 
simulation). 

The comparison is done in terms of several different metrics, including 
realized detection costs c(XQ..,.{t)), distribution of detection times r, and fre¬ 
quency of false alarms, represented by d{Xr) = 1 — Pr in our setup. As 
expected, the Optimal strategy with detection time r* that directly opti¬ 
mizes the cost-benefit in the full model performs best. The corresponding 
expected costs are F(xo) ^ 6.53, with average detection time E[t*] ~ 8.86. 
It outperforms the Threshold-P strategy by about 7% in terms of reducing 
detection costs, and the Threshold-t strategy by about 9%. These are non¬ 
trivial cost savings which highlight the benefit of information fusion. Table 
also shows that the 2-D LP approximation performs well in this example, 
generating very similar expected costs. At least for this case study, detec¬ 
tion happens early enough that the branching process approximation of the 
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outbreak works fine. 

Recall that our model is stochastic and generates adaptive detection 
strategy. Hence the detection time t* is a random variable. As shown in 
Table the corresponding standard deviation StDev{T*) ~ 2.6 is sub¬ 
stantial. This illustrates the sub-optimality of the Threshold-t strategy 
that stops at a fixed t with StDev{t) = 0 trivially. Not surprisingly, the 
ability to delay or speed up outbreak announcements based on latest data 
are crucial for optimizing policy making. We also note that compared to 
the Threshold-P strategy, the Optimal strategy tends to announce later, 
E[t*] ~ 8.86 > 7.88 ~ this is also confirmed by the respective 

histograms of t* and in Figure However, we emphasize that the 

detection rules do not have a clear ordering. In other words, the random 
variables r*, etc., cannot be directly compared. 

A complementary metric of detection quality is provided by the prob¬ 
ability of false alarms, PFA := £'[1 — Pr\. For the optimal strategy we 
find that PFA* = 8.2%. In contrast, for Threshold-P strategy, we have 
ppj^Thr-P _ Note that because we use a discrete-time model, 

at time of detection Pr will strictly exceed the threshold P = 0.8, hence 
< 1 — P. The histograms of Pr are shown in Figurej^and confirm 
the qualitative difference among the detection strategies. The Threshold-P 
strategy only stops once Pt > P, so that Pr has support on roughly [0.8, 0.9]. 
In contrast, the adaptive Optimal (and LP) strategies, have a much wider 
range for Pr- In particular, sometimes epidemics are announced even before 
Pt hits the level 0.8. 

To further quantify the improvement provided by the Optimal detection 
rule, Figure gives a scenario-by-scenario comparison of relative realized 
detection costs. Note that in hindsight, t* may sometimes perform worse 
that or even Figurej^plots the histogram of the difference in 

costs for each trajectory Xg.^, n = 1 ,..., 1000, namely c(xo:t*)) c(xg..,_T/ir-p), 
and c{xQ.rThr-t). We find that the costs computed with Optimal/LP strate¬ 
gies are smaller than costs computed with Threshold strategies for more 
than 80% of the trajectories. 

To sum up, we observe material improvement from using Optimal de¬ 
tection rule in this case study. Moreover, the obtained detection rule is 
substantially different from the thresholding protocol. On the one hand, 
the adaptive detection time t* exhibits a wide spread and is highly non¬ 
constant across trajectories. On the other hand, the posterior probability 
of false alarms Pr* is also strongly variable. As a result, the average fre¬ 
quency of false alarms is drastically lowered relative to Threshold-P strategy. 
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Difference in costs 

Figure 9: Relative realized detection costs across different strategies. 
The histogram shows the distribution of the difference in costs along the 
1000 simulated trajectories, namely c(xo;r*) — c{'x.Q.^Thr-p), and c(xo:r*) — 

c{^Q.^Thr — t ) . 


reducing overall expected costs. 

4-2. Effect of Detection Cost Parameters 


Cfa 

Mean 

T* 

StDev. 

Cost 

Mean StDev. 

PFA = E[1 - Pr*] 

10 

6.84 

1.62 

5.32 

0.99 

21.4% 

20 

8.87 

2.60 

6.54 

1.71 

8.3% 

30 

9.61 

2.79 

7.21 

2.22 

5.3% 


Table 3: Summary statistics of the Optimal detection strategy r* for differ¬ 
ent false alarm penalties Cfa- 


The main parameter in our quickest detection setup is the ratio of the 
cost of false alarms and the cost of detection delay, C'pA/C'Deiay A high 
ratio penalizes premature announcements and requires more care in the 
assessment of the potential outbreak in Pool 2. A low ratio invites more 
aggressive actions. To better understand the role of this ratio, in Figure [l0| 
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Number of Infecteds in Pop. 1 


Figure 10: Boundaries of detection maps constructed based on differ¬ 

ent penalties for false alarm, Cfa- 


we show several detection boundaries corresponding to varying Cfa, 

while Coeiay = 1 is kept fixed. As expected, a lower Cfa enlarges the 
Announce set ©. In particular, the boundary d& shifts down and to the 
right. As a result, starting from a fixed location the stopping set 

& will be reached sooner, so that r decreases (in the sense of stochastic 
dominance for the corresponding random variables). This is confirmed in 
Tablej^that reports statistics for r* and various Cfa- We hnd that E[t*] = 
8.86 when Cfa = 20, but is only E[t*] = 6.84 for Cfa = 10. Simultaneously, 
the frequency of premature announcements PEA will increase. The precise 
relationship is however nonlinear. Lowering Cfa from 20 to 10, the PFA 
rises dramatically to about 21% from 8%. Conversely, raising Cfa to 30 
only reduces PFA to 5.3%. A common approach in the decision literature 
is to select a priori a desired level of PFA (say PEA = 10%) and then 
numerically solve the inverse problem to obtain the corresponding Cfa and 
hence the corresponding detection rule &. 


5. Numerical Implementation 


, defined recursively in 
(10) we use approximate dynamic programming techniques. In particular. 


To hnd the detection maps 6 t for t = 1,2, 
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we rely on the Regression Monte Carlo approach miE] to approximate the 
conditional expectation map over x G Al. 


5.1. Regression Monte Carlo 

For the remainder of this section the auxiliary “time” variable t is fixed 
and the goal is to approximate the conditional expectation q{t, x) := Fl[c(XQ.^(t) |Xo 
x] in 0- Recall that at step t, detection rules are restricted to satisfy 
< t. The Regression Monte Carlo technique approximates q{t, •) by a 
predicted surrogate value q{t, •) which is based on a statistical regression 
framework. 

The surrogate prediction is built using data simulated from the specified 
model. To do so, a design Z := {xq, n = 1,... ,N} of N locations is first 
generated. Next, we generate the corresponding scenarios {Xg.^} with the 
initial value Xg = Xg, one scenario for each initial location. Define 

r” := min{s > 1 : X” G (19) 

which leads to path-wise waiting costs g” := c(XQ..^n) using formula Q on 
the n-th scenario. The aggregate dataset is 


Z = {(x^,g"),n = l,...,iV}. (20) 

The construction of q{t, •) then involves response surface modeling, i.e. de¬ 
termining the relationship between the initial condition x and the mean of 
the sampled Q\x = c(Xo:rt)- Statistically, we start with 

Q\x = q{t,x) + e, (21) 


where q{t, •) is the true response surface, Q = c(Xg..,.(t)) are random scenario- 
based costs, and e are mean-zero residuals with variance cr^ arising from 
Monte Carlo simulations. Empirically, (21) translates into regressing {g”} 
on {xg}, n = 1,..., A^; this step is discussed in section 5.2, After determin¬ 
ing q, and using (10) the estimated detection rule &t is 


&t := {x : q{t,x) - d{x) > 0} . 


( 22 ) 


The above provides a recipe to obtain an (approximate) &t using the 
collection of detection rules Iterating over t, yields the sequence of 

detection maps &t for t = 1,2,.... We recall that as t —)■ oo, we expect &t 
to stabilize and tend to a time-invariant detection map. Such convergence is 
illustrated in Figure [ll| where we trace the boundaries d&t for t = 1 ,..., 20. 
Convergence takes hold after about 15 iterations and suggests that ©20 — ©; 
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this is what we used for Figures hflO where the boundary of 620 was taken 
as the final output of the Algorithm. 


The detection rule (19) is time-dependent since it utilizes a new &t-s 
at each stage s. Model predictive control simplifies this feature with a 
time-invariant rule that simply utilizes &t-i (that we relabel as 
for typographical distinction). Indeed, as Figure 11 shows, the early maps 
Si, 02 ,..., are not as accurate as &t-i for t large, so it makes sense to com¬ 
pletely “forget” them and rely just on the last iteration step. Accordingly, 
we implement a blend of 0 and ([^ by first using ( [T^ over t = 1,2,... ,t* 
and then switching to a receding-horizon rule 


^MPC ._ > 1 : Ts G Si_i}, t = t*,t* + 1,.... 


(23) 


The above MFC iterations are terminated once q{t,x.) and g(t-|-l,x) do not 
change much, namely ||g(t, •) — q{t - 1 - 1 ) •)IU°° < Tol for a specified tolerance 
level Tol. 



Figure 11: Convergence of the detection boundaries d&t over t = 1 to t = 20 
for the 2-D LP detection rule from Section 


5.2. Regression Model 

Because we have limited a priori knowledge about the structure of the 
detection rule, it is preferable to work with a nonparametric regression ar¬ 
chitecture for g(t, x). (For example a linear regression model for q would 
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imply that S in (22) is defined through linear constraints, i.e. forms a sim¬ 
plex m X.) In addition, nonparametric regression is typically more robust 
for dealing with the non-Gaussian residuals e that arise in our model. 

There are numerous nonparametric regression frameworks that can be 
used, including splines, Gaussian processes, or generalized additive mod¬ 
els; see e.g. the classic monograph [35]. Note that even though x i—>■ ^(x) 
is continuous, some discontinuous response surfaces might also be helpful, 
such as random forests or dynamic trees |9]. In the present article we take 
up a variant of local linear regression, known as Loess. Loess fits weighted 
linear regression models to localized subsets of data, determined using a ker¬ 
nel function, specifically a fe-nearest-neighbor algorithm |3S|- Gompared to 
classical linear models, Loess better handles outliers and heteroscedasticity, 
and also does not make assumptions about the global shape of the response 
surface. 

The Loess response model is of the form 


9(Lx) = ^/3i(x)fij(x) 
i=l 


(24) 


where Bi(-) is the set of r pre-specified basis functions and /3i are estimated 
regression coefficients at x. Given input matrix X and matching response 
vector Q, /3 is fitted using local least-squares minimization 

/3(x) := arg min Kx{^,X){Q - B{X)'^0f, (25) 

/3eK’’ 


where Kx(^x.,X) is the weighting kernel. The idea behind the kernel is to 
base the predicted q{t, x) on the samples in the neighborhood of x, weighted 
by their distance from x |45[ Sec. 2.8.2]. The size of the neighborhood is 
controlled by the smoothing parameter A. If A < 1, only a proportion A of 
the samples will be used in fitting. The smaller A, the more “wiggly” the fit 
q{t, •) is going to be since fewer samples are used in computing /3(x). Loess 
can be viewed as a special kernel regression method, with the prediction 
being a weighted average of the responses q"^: q{t,x) = for the 

equivalent kernel Z(-). In our numerical examples, we use the implementation 
of Loess provided in the R by the built-in package stats [37|, which uses a tri- 
cubic kernel and linear, first-order basis functions; the smoothing parameter 
was A = 0.4. 


5.3. Experimental Design 

The aim of the response surface is to maximize the accuracy of Sf. 
This is equivalent to maximizing model fidelity along the boundary of the 
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detection map. Statistically, for a localized response surface, accuracy is pri¬ 
marily driven by the local density of the input data that is specified by the 
experimental design Z. Hence, to maximize our confidence regarding the 


boundary of ©t in (22), we generate appropriate, adaptively chosen exper¬ 


imental designs Z. This is achieved using the Sequential RMC framework 
introduced by [9]. SRMC uses tools from active learning/Bayesian optimiza¬ 
tion to gradually grow the design Z so as to zoom-in to the boundary of 
This is done by first quantifying the accuracy of the existing response 
surface, and then adding new design sites so as to maximize information 
gain. See for details. The SRMC approach is illustrated in Figure 

where the adaptively generated experimental design Z (of size 2000 in 
the hgure) is highly concentrated around the detection boundary d&. This 
targeted sampling of outbreak scenarios allows for more efficient estimation, 
in particular lowering the local standard errors 'O(x) along d&t, cf. the right 
panel of Figure]^ 

In ( |22[ ) the boundary of &t corresponds to the regions of X where the 
cost difference between immediate detection and waiting is zero. Hence, we 
aim to have more design points in regions where {g(t,x) — d{x.) ~ 0}. To 
this end, we define the “posterior” measure of response surface accuracy via 


p(x) := 


-\q{t,x) - d(x)| 




X 


(26) 


where is the standard normal cdf and the predictive variance v measures 
the standard error of the surrogate prediction. 


u(x) = (T^(x 


X 


(27) 


with (T^(x) the estimated variance of e around x in (21), see |45L Sec 6.1.2]. 


The motivation for (26) is that p(x) mimics the Bayesian posterior prob¬ 
ability of estimating the wrong sign (conditional on the samples in Z) of 
q{t, x) — d(x), assuming that the posterior distribution is Gaussian with the 
empirical mean q{t,x) and variance ^(x). 

The defined metric p{-) serves as a guide to augment new design loca¬ 
tions. Namely, it defines an acquisition function w(x) for greedily growing 
Z, similar to active learning methods j49| . The acquisition function is high¬ 
est in the regions where p(x) is close to 0.5 which correspond to d&f Our 
main choice is 


y^mm(x) _ [p(x), 1 — p(x)] . 


(28) 
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Alternatives include the Gini weights = p(x) (1 — p{x)) and En- 

tropic weights = —p(x) logp(x) — (1 — p(x)) log(l — p(x)). 

To speed up the response surface modeling, which requires refitting of 
q{t, •) multiple times, we used batch steps, incrementally working with de¬ 
signs of size N = Nq, Nq+N',..., At each sequential design iter¬ 
ation, an additional N' design points are added to existing 

Those are sampled multinomially in proportion to the acquisition function 
w{-) from a candidate set Aignite. Both the initial design and the 

candidate sets Xfinite are generated using Latin hypercube sampling (LHS) 
of size D from X. The overall procedure, summarized in Algorithm A..2\ 
finally refits at each iteration the Loess model for q (and hence 6^), grows 
the experimental design = Z^^'^ U {^o}n=N+i recomputes the 

acquisition function (28). As the design size gets larger, we expect that the 
implied empirical estimate gets closer to the true d&t- 


Remark 3. One can apply standard, nonsequential RMC by skipping the 
inner while loop (steps 7-15) in Algorithm A. A This reduces to building a re- 
sponse model on a pre-specified (possibly randomized) design Z := {xgj^^j^, 
keeping all other steps as is. 


For the detection map in Figurein Section]^ we used an initial design 
of A"o = 200, which was grown over 10 iterations with N' = 200 to a final 
design of = 2000. The acquisition function was tc™™ and the candidate 
sets Aifinite of size D = 2500 were generated with LHS. Since detection 
happens while is still relatively small, we restricted the response surface 
regression domain to G {0,1,..., 400}, G {1000,..., 2000}. Lastly 
we note that the method is still computationally intensive, with the bulk of 
the effort spent on generating T ■ scenarios of X; running times (on a 
8-core 2.27GHz machine with 12GB of RAM) were about 20 minutes. 


6. Discussion 

We have presented a framework for optimal detection of epidemics in a 
coupled meta-population model. Our approach explicitly takes into account 
cost-benefit considerations regarding announcement of an epidemic, as well 
as spatial dependence across susceptible pools. Given the information about 
two populations and characteristics of the infection, our algorithm produces 
the full detection map which can then be used repeatedly. We demonstrate 
that information about the epidemic in one pool can be used to lower the 
detection costs in another pool, realizing savings compared to traditional 
threshold-type detection methods. 
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Since our dynamic optimization approach is entirely simulation-based 
it is unusually flexible. Indeed, the precise underlying epidemic model of 
X is not crucial, since Algorithm |A.2| only requires the ability to generate 
its trajectories. In fact, the computational complexity of our algorithms 
is tied not to the dynamics of X but to its dimensionality. In Section 
we had dim(X) = 3; based on our experience with RMC in [301 E]) 
present approach can straightforwardly handle up to 6-8 dimensions. In 
high dimensions, extra care must be applied for generating the experimental 
designs Z since the concept of neighborhoods underlying nonparametric 
local regression breaks down. For example. Loess regression performs poorly 
if the dimension of the data is higher than 4-5. 

The presented SIR framework gives a basic mechanistic description of 
disease progression that is obviously not very realistic. More sophisticated 
versions might allow for further compartments (such as Exposed or Dis¬ 
eased individuals), age stratification, and heterogeneity among the meta¬ 
populations. One could also include further transitions beyond (13), such 
as immigration 0 —)■ immunity lapses —)• or vaccination 

—7- Introducing immigration would allow for endogenous epi¬ 

demic in Pool 2, removing the assumption that outbreaks always start in 
Pool 1 and then spread to Pool 2. The constant transition rates used can be 
replaced with seasonal patterns, stochastic fluctuations |50j . or hierarchical 
Markov structures [2]. 

Alternatively, one can also imagine more sophisticated models for the 
outbreak pseudo-posterior Pt - recall that the proposed one was largely for 
convenience than any realism. For example, the Gaussian noise 6t in the 
dynamics of Pt that was used in the case study may be better modeled via 
a Beta distribution (which arises naturally as conjugate to the Poisson in¬ 
crements of the fully-observed stochastic SIR model |2]). Overall, the key 
requirement is the Markov structure which makes it possible to use regres¬ 
sion against X to describe the detection rule. The Markovianity requirement 
can be partly relaxed if one is willing to accept approximately-optimal so¬ 
lutions. Indeed, one can always project the optimal detection rule into the 
smaller space of rules that only depend on some subset X'; in other words 
restricting the detection map to only take into account some of the state- 
space dimensions. This idea was already discussed in Section [^ where we 
described the sub-optimal LP strategy. 

Second, one may modify the cost structures (®-@ to better capture 
the desired detection goals. The presented costs were motivated by their 
classical analogues in sequential change-point detection, but might not be 
the most appropriate for public health contexts. For example, there is some 
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leeway in what constitutes an outbreak. In Q, the threshold was zero, 
i.e. even a single infected individual in Pool 2 was reportable. One can use 
thresholds I other than zero, so that an outbreak is reportable only when 
II ^ > I, otherwise an announcement is treated as premature. Similarly, 
the waiting cost in ® was constant; it may be more realistic to make it 
( 2 ) ^ 

proportional to ^, which would correspond to fixed costs per infected. 

For such more general formulations, the costs d{X) and c(X) would no 
longer be functions of Pt, and one would need to work with the full posterior 
distribution tt* of \Qt- The RMC framework could still be usable, namely 
we may use particle filtering m to obtain vr^ along a simulated trajectory 
of the underlying epidemic model. Certainly, particle filtering can become 
computationally expensive, making efficient inference essential. We refer 
to mm for some recent implementation strategies in this direction that 
specifically target epidemic models. The integrated sequential inference + 
optimization model would then allow to treat a partially observed version 
of a iiT-pool SIR model of (13), and ultimately a larger-scale setup such as 
influenza surveillance across all 50 states, cf. Figure 
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Appendix A. Algorithms 


Algorithm A.l Path and Cost Generation 

Require: 6o:i-i 

1: for n = 1,..., A do 
2: S ^ 1 

3: while s < t do 

4: Simulate the next state x” ~ pi(-|x”_^) 

5: if Xg G &t-s then Break 

6: end if 

7: S S + 1 

8: end while 

9: rf s 

10: Compute (?"■ = c(xq.^„) using formula Q 

11: end for 

12: return {(x[J, 
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Algorithm A.2 Sequential Regression Monte Carlo 


Require: C'FA,CDeiay, Aq, A', D 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


6o ^ A 

for t = 1, 2,... do 

Generate experimental design {xq, 
Compute scenario costs = c{Xq.^ 

rithm A.l and ©o i-i 


n = 1,..., Ao} 

n) for n = 1,..., Aq using Algo- 


Regress {(?"'} on {xq}, n = 1, 
Initialize A •(— Nq 

while N < do 


, Aq using Loess (24) 


Generate Agnite of size D using Latin Hypercube Sampling on X 
Compute the acquisition weights w{x) Vx G Agnite via (28) and 


(26) 


Sample {xo}^+^_^^ from Agnite using weights w{x) 

Simulate the costs g”, n = A+1,...,A+A' using Algorithm |A.1| 

Update the Loess regression model (24) using the latest Z 
N ^ N + N' 

end while 

©t {x G A : q{t, x) — d(x) > 0}, cf. (22) 

end for 
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